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We examine an extraction model for metamaterials, not previously reported, that gives precise, 
quantitative and causal representation of S-parameter data over a broad frequency range, up to 
frequencies where the free space wavelength is only a modest factor larger than the unit cell dimen¬ 
sion. The model is comprised of superposed, slab-shaped response regions of finite thickness, one 
for each observed resonance. The resonance dispersion is Lorentzian and thus strictly causal. This 
new model is compared with previous models for correctness likelihood, including an appropriate 
Occams factor for each fit parameter. We find that this new model is by far the most likely to 
be correct in a Bayesian analysis of model fits to S-parameter simulation data for several classic 
metamaterial unit cells. 


I. INTRODUCTION 

The field of metamaterials has promised a dramatically 
expanded range of material properties to the electromag¬ 
netic designer. However, the compelling performance 
gains that could be realized in many devices, are tem¬ 
pered by two problems. One is that the desirable range 
and tunability of the real parts of the permittivity and 
permeability come also with undesirable qualities: loss 
and spatial dispersion. Even quantifying these undesir¬ 
able qualities is a challenge. Not only is there no simple 
standard metric for quantifying spatial dispersion, but its 
presence makes dubious the practice of quantifying loss 
with material property imaginary parts. 

The second problem is that we lack robust algorithms 
for assessing a given metamaterial design without ad-hoc 
human intervention. Ideally, our algorithm would pro¬ 
vide simple and physically meaningful, quantitative de¬ 
scriptions of the effective behavior. Without such algo¬ 
rithms one cannot hope to exploit the large scale compu¬ 
tational resources and optimization strategies that would 
lead to superior designs. 

The premise of the current work is that extracting 
metamaterial properties from simulation data into the 
best model will provide progress with both of these prob¬ 
lems. For the purpose of this article, the best model is 
the simplest one that provides accurate quantitative rep¬ 
resentation of the simulated behavior. Signihcant spatial 
dispersion seems to be unavoidable with metamaterials 
that provide their unique (and sometimes extreme) prop¬ 
erties with practical unit cell dimensions. Several authors 
have suggested that a spatially dispersive model provides 
more physical insight than a spatially local on^ilt^. In this 
work we describe several models that incorporate spatial 
dispersion through unit cell inhomogeneity. We analyze 
these inhomogeneous models - including some not previ¬ 
ously discussed for metamaterials - in an objective sta¬ 
tistical analysis to identify a preferred model for several 
typical unit cell designs. 

All the extraction models presented here are comprised 
of homogeneous, isotropic, layers, for which the reflection 
and transmission coefficients of normally incident plane 


waves can be analytically calculated. A metamaterial 
unit cell is simulated and its normal incidence reflection 
and transmission coefficients computed. Model parame¬ 
ters are found by fitting reflection and transmission co¬ 
efficients from the model to those of the simulation data. 
These model parameters include the thicknesses of su¬ 
perposed, slab-shaped response regions, and their corre¬ 
sponding Lorentzian susceptibilities. The fits are per¬ 
formed over the entire frequency range simultaneouslj^^. 

Of note, we do not incorporate spatial dispersion or 
magneto-electric coupling into the slabs response itself, 
though in a homogenized sense, the layered structures 
can represent behavior that appears as such. If the 
models perform well, then adding such complexity is 
unnecessar 

The simulations presented here have been performed 
with the FDTD solver of CST Microwave studio, but 
other solvers should provide similar results. 


II. THE FOUR MODELS 

In describing the extraction models we will analyze, 
we will define the concepts of response slabs and mate¬ 
rial layers in a specific way. Geometrically, both con¬ 
structs are assumed to be infinite in extent in the direc¬ 
tions transverse to the wave propagation direction. Re¬ 
sponse slabs have a thickness and a spatially uniform 
susceptibility (either electric or magnetic) arising from 
a single resonant mode. Since the total response at any 
given point in space is the sum of the responses due to all 
resonant modes that extend to that point, we represent 
the total response by a superposition of response slabs. 
This superposition of slabs creates a number of distinct 
uniform material layers with uniform material proper¬ 
ties. In a given layer, these properties are just the sum 
of the susceptibilities of the slabs that overlap that layer. 
For example, in FIG. 0a) , the unit cell has two over¬ 
lapping slabs, but five material layers. Two outer layers 
have vacuum properties (zero electric and magnetic sus¬ 
ceptibility). Moving inward, there are two layers with 
the material properties give by only the susceptibility of 
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FIG. 1. The four models: (a) multi-thickness, (b) single¬ 
thickness, (c) thin sheet, and (d) homogeneous. Two slab¬ 
shaped response regions (red and green) are shown, but an 
arbitrary number are allowed, one for each resonance. (Slabs 
are of infinite transverse extent.) 


the green slab, and a single, central layer with material 
properties equal to the sum of the susceptibilities of the 
green and red slab. For this work we assume that the re¬ 
sponse slabs are centered in a cubic unit-cell of dimension 
d, and each slab, i, with thickness, si, is associated with 
an electric or magnetic response whose unit-cell-averaged 
susceptibility, Xi, is described by a Lorentzian dispersion, 

Xz (/) = ( 1 ) 

The Lorentzian dispersion has the usual parameters of 
static response, resonance frequency, fi, and loss, Si. 
In this function, we do not include the high-frequency 
response (x“) which arises from resonances above the 
frequency range of interest. Instead, we include any re¬ 
sponse due to above-range resonances as a separate slab, 
with its own thickness. The susceptibility of such a slab is 
the constant, Xi^ fh® high resonance-frequency (i.e. low 
frequency) limit of 0 . This treatment is motivated by 
the fact that above-range resonances won’t necessarily 
have the same mode shape as any of the in-range reso¬ 
nances. It should be noted that we have allowed a con¬ 
stant magnetic susceptibility to be negative in all three 
of the examples anal^ed in this article. As pointed out 
by Wood and PendrjffI, this can be allowed in causal me¬ 
dia. Below we highlight the one case where causality has 
been compromised in a specific fitted model. 

The unit-cell-averaged susceptibility is defined in re¬ 
lation to the local, slab susceptibility by equating the 
polarizabilities associated with the appropriate volumes. 

< = ^ = £oX^s^d^ = ^oXid^ (2a) 

aT='^=xTs.d^^xTdP (2b) 

For either electric or magnetic response, we find the lo¬ 
cal, slab susceptibility is related to the unit-cell-averaged 
susceptibility through the thickness ratio. 

Xi=Xi— (3) 

Si 

The reason to associate the Lorentzian parameters with 
the unit-cell-averaged susceptibilities instead of the lo¬ 
cal, slab susceptibilities is to decouple the thickness and 
Lorentzian parameters. For example, in this scheme, the 
limit of zero thickness (with finite polarizability) can be 
treated without divergent Lorentzian parameters. Also, 
we prefer susceptibilities to permittivities and permeabil¬ 
ities since their superposed response can be calculated by 
direct summation. The four models we consider are de¬ 
scribed below. We begin with the most general, with 
each Lorentzian slab allowed its own unique thickness. 
The other models are special cases that follow from the 
general one: a single thickness for all slabs, all slabs of 
zero thickness, and all slabs of unit cell thickness, (i.e. 
the homogeneous model). 
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A. Slabs of multiple thicknesses 


The superposition of the susceptibilities of N overlap¬ 
ping slabs centered in the unit cell, creates up to 2N + 1 
distinct, contiguous, uniform material layers, including 
two vacuum layers if the thickest slab is thinner than the 
unit cell dimension. The refractive index and impedance 
of these layers is given by 
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where the index, i, refers to the slabs, and the index, j, 
refers to the layers. The notation, i S j, means sum over 
all of the slabs, i, that overlap in the layer, j. We assume 
these contiguous layers are ordered from front-to-back. 
The front of the unit cell is the reference plane for the 
reflection coefficient, r,and the transmission coefficient, 
t' , refers to complex plane-wave amplitude at the back 
of the unit relative to the fronf following the notation 
of Smith. We can find these coefficients using transfer 
matrix methods. In particular the transfer matrix of the 
individual layers is given by 


T,= 


Uj 


aj = cos irijklj) 
13j = sin {ujklj) 




iPjCj aj -b i^jC 




(5a) 

(5b) 

(5c) 

(5d) 


where Ij is the thickness of a layer and k is the free-space 
wave-vector. The transfer matrix of the combined set 
of layers just the matrix product of the individual- layer 
transfer matrices. 


T = nT, (6) 

j 


and the reflection and transmission coefficients are found 
simply from the elements of the combined transfer ma¬ 
trix. 


r= ^ 
Tn 

ill 


(7a) 

(7b) 


Unless there is just one slab, it is clear that equations 
Q cannot be inverted to find the material properties of 
the slabs as a function of the reflection and transmis¬ 
sion coefficients independently at each frequency point, 
since the number of unknowns exceeds the number of 


equations. Instead, to find the slab material properties 
from simulation derived reflection and transmission co¬ 
efficients, we perform a least-squares fit over a range of 
frequencies, employing a distinct Lorentzian dispersion 
for each slatP. There are two or four fit parameters per 
slab depending on whether the slab represents an above¬ 
range resonance or an in-range resonance, respectively. 
The full set of fit parameters is 


si,Xi 

or 

si,Xi,/i,^i 

S2,X2 

or 

S2,X2:/2,^2 

sn,X%[ 

or 

SN,XN’ 


If the number of slabs is appropriate to the number of 
resonances in the range of frequencies used, there is suffi¬ 
cient information in the reflection and transmission data 
to precisely determine all of these parameters. 


B. Slabs of a single finite thickness 


The single slab is the same as the multi-slab, except all 
of the Lorentzian response functions are associated with 
a single thickness, s. The refractive index and impedance 
for the single material layer that is coincident with this 
slab is a simplification of equations Q, 



z = 


i + fEitr 

i 

i + !eE' 


(9a) 


(9b) 


In addition to this slab, the unit-cell is comprised of two 
surrounding, free-space layers of thickness (d — s)/2. It 
is convenient to define reflection and transmission coef¬ 
ficients, rs and t(,, that are referenced (de-embedded) to 
the front and back of the slab. 


t' = 


(10a) 

(10b) 


These new coefficients can be found from a single transfer 
matrix, which results in the frequently used expression^. 


1 . , ^ 1 ■ 
— = cos [nks) + -J 
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The equations (11) can, of course, be inverted to yield 
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(12a) 

(12b) 
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FIG. 2. Tuning the model slab thickness, s, to a physically 
preferred value that minimizes non-causal response. The unit 
cell is a split ring resonator (microwave cloak, design 1, in¬ 
set). The frequency range shown includes the lowest electric 
resonance (above the lowest magnetic resonance). The blue, 
red, yellow and green curves are for s = 3.33 (full unit cell), 
2.54, 1.74 and 0.95 mm, respectively. At the preferred thick¬ 
ness (green curve) the permittivity approaches the Lorentzian 
form and the anti-resonance response of the permeability has 
been eliminated. Note that for thicknesses thinner than the 
preferred value, non-physical features return to the response 
(not shown). 


SO that the slab material properties can be found from 
simulation derived reflection and transmission coeffi¬ 
cients, independently at each frequency. We can also 
And the unit-cell-averaged susceptibilities. 






= ( 1 ) ^ 


(13a) 

(13b) 


The slab material properties, , as well as unit-cell- 
averaged susceptibilities, (13), are dependent on the 
value chosen for the thickness, s. One can choose this 
thickness to minimize spatial-dispersion induced arti¬ 
facts. An example is shown in Figj^ 


C. Thin sheet: Slabs of infinitesimal thickness 

Consider the case of a super-position of slabs, all of 
infinitesimal thickness, but where we assume the unit-cell 


averaged susceptibilities are flnit^l^. To take the thin- 
slab limit of equations ( [l0| , 0 and ( |12[ ) , we require the 
following limits, which we And using equations (§. 


lim ns = 


lim z = 

s—>0 
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X/ Xrr, 
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= ZO 


' Xmi d = nod (14a) 


(14b) 


As seen from ([^, the local susceptibilities are diver¬ 
gent, so that the unit constant in (|^ can been ne¬ 
glected by comparison. The limits are, however, finite. 
On the right-hand-side we define special index, uq, and 
impedance, Zq, for this zero-thickness layer. We find the 
relationship between the thin slab coefficients tq and tg 
and the de-embedded coefficients r and t' by taking the 
same limit of equation (101 


r = e-^“ro 
t' = 


(15a) 

(15b) 


Additionally, taking the limit of equations (11) and (12) 
we find 


^ = cos {nokd) + ]-j 
to 2 


zo 


— ] sin (nokd) (16a) 
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(16b) 


(17a) 


(17b) 


We can also include infinitesimal slabs in a multi¬ 
thickness slab. All of the overlapping infinitesimal slabs 
will contribute a single transfer matrix, T^-. This transfer 
matrix is found by taking the Ij —> 0 limit of equations 
([5]), just as we took the limit of equations ^ and ( [T^ . 
Any finite slabs that overlap the infinitesimal slabs do not 
contribute and can be neglected for this transfer matrix. 


D. Homogenous: Slabs all of thickness d 


The standard model for material property extraction 
is the uniform or homogeneous medium. The primary 
advantage of this model is simplicity. Engineers already 
know how to design and analyze devices that incorpo¬ 
rate uniform or slowly varying media. Like the single¬ 
thickness and thin sheet models, the homogeneous model 
can be inverted, and the inverse formulas applied inde¬ 
pendently at each frequency point. Equations (9p2) can 
be applied by setting s = d. 
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FIG. 3. The three unit cells analyzed: (a) the microwave in¬ 
visibility cloak, cylinder 1 design, (b) a two-dimensional neg¬ 
ative index medium design, (c) the original ELC resonator 
design. 

Obviously, these material properties will precisely rep¬ 
resent, with the use of the Fresnel formulas ( [IT| ), the 
reflection and transmission behavior for a sample of the 
same thickness from which they were extracted (though 
not very usefully). They may also precisely predict 
the reflection and transmission behavior for an arbitrary 
sample thickness, if the inter-cell, electromagnetic cou¬ 
pling is restricted to dipole interactions. Also, if a suf¬ 
ficient number of unit cells are used in the propagation 
direction, the material properties will converge to those 
of a bulk medium, regardless of the coupling behavior. 
The extracted bulk properties will apply to all sufficiently 
thick slabs. 

Unfortunately, these material properties will usually 
provide limited physical insight, because the model is of¬ 
ten too simple to represent the metamaterial in ques¬ 
tion. Wedging the reflection and transmission behavior 
into this over-simplified model leads to material prop¬ 
erties that violate physical principles. Specifically, the 
practical, moderate values of free-space phase advance 
across the unit cell (i.e. Ao/d), together with the extreme 
polarization found near a resonance, lead to significant 
field variation across the unit cell, which is not compat¬ 
ible with an effective, uniform (spatial dispersion free) 
medium. The resulting material properties are then a 
poor fit to the simple Lorentzian line sh apes o ne expects 
from an isolated, second-order resonanc c^^ ^ ^^ l 


III. EXAMPLES 

We demonstrate the fitting of the four models to sim¬ 
ulation derived reflection and transmission coefficients 
(S-parameters) for two well known unit cells: the mi¬ 
crowave invisibility cloak (Figj^)P^ and the ELC res¬ 
onator (Fig(^)P. Another unit cell, one designed for 
two-dimensional, isotropic negative index is also included 
(Figj^). This latter design provides a somewhat more 
complex response with three overlapping resonances. 


A. Microwave Cloak Unit Cell 

Fitting of all the models provides good constraint of 
the fit parameters, as seen from the parameter uncertain¬ 


ties in Table [Tj However, only the multi-thickness model 
provides a quantitatively precise fit. In plots covering the 
entire fitted frequency range (such as Fig. |^a) and (b)), 
simulation and model curves would be indistinguishable. 
The small deviations between model and data can only 
be observed in zoomed-in sections of the plot (Fig. Qc)) 
or in plots of the residuals (Fig. |^d)). The chi-squared 
per degree-of-freedom and probability measures of model 
appropriateness give quantitative support to choosing the 
multi-thickness model. (These measures are normalized 
to the multi-thickness model fit, as described in section 
IV.) The probabilities for the other models are so small 
that they are essentially zero. 

One subtlety arose in obtaining the best fits for the 
multi-thickness model. As can be seen in Table |Tj three 
of the four thicknesses are given as approximately zero. 
Here the results for a zero thickness slab (equations ( |l4| )- 
0 ) were not used. Instead, a different small but finite 
thickness was used for each slab, with Si < S 3 < S 4 . 
The quality of the fit was independent of these thickness 
for any so-ordered set of values less than about 0 . 01 mm. 
Good, but less impressive, fits could be obtained using a 
zero-thickness slab in the multi-thickness model. 


B. 2D Isotropic Unit Cell 

As with the cloak unit cell, here only the multi¬ 
thickness model provides a quantitatively precise fit as 
seen from the chi-squared per degree-of-freedom and 
probability measures of model appropriateness in Ta¬ 
ble [H] Particular to this unit cell, the homogeneous 
model obtains the same quality of fit while omitting 
the second electric resonance included in the other mod¬ 
els. This is more an indication of a poor fit, than ev¬ 
idence for the absence of the resonance. Though the 
multi-thickness model includes fourteen parameters—a 
substantial number—the complexity of the four, real¬ 
valued data sets (i.e. the real and imaginary parts of 
the reflection and transmission curves) is such that only 
a physically-motivated, and correct model has any chance 
of providing a good fit. This is born out in the proba¬ 
bility measure, which includes an Occams penalty factor 
for each parameter (as describe in section IV). Here the 
zero-valued thicknesses in the multi-thickness model do 
use the limiting results of the zero thickness slab, equa¬ 
tions 

C. ELC Unit Cell 

For this unit cell, Fig|^c), the results are much the 
same as the other two, in terms of the quality of the 
model fits. However, we point out two problems that 
occurred with the multi-thickness model fit, as shown in 
Table |HI[ One is the that thickness of the constant mag¬ 
netic susceptibility slab, S 3 , exceeds the dimension of the 
unit cell, d, and all of its material components. This 
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FIG. 4. Microwave cloak unit cell. From simulation, the real and imaginary part of the (a) reflection coefficient and (b) 
transmission coefficient, (c) Best global fits of the four models to the simulated data (black dots). The narrow range of 
frequencies displayed is indicated by the dashed box in (a). The models are shown in: green (multi-thickness), yellow (single 
thickness), red (thin sheet) and blue (homogeneous), (d) The combined residuals, A, as defined in equation (26l, with the 
same color scheme. 


TABLE I. Microwave cloak unit cell. 



multi-thickness 

single thickness 

thin sheet 


homogeneous 

xVdof 

1 

80 

200 


1000 

probability 

1 

~ 0 

~ 0 


~ 0 

parameters 

9 

9 

8 


8 


electric parameters 


Si 

X ? 

h 

Si 

~ 0 

0.31754 

21.6825 

0.0506 


± 0.00004 

± 0.0002 

± 0.0002 


1.2091 

0.36051 

21.4801 

0.0587 


± 0.0004 

± 0.00006 

± 0.0002 

± 0.0002 


0 

0.26109 

21.9129 

0.1911 


± 0.00007 

± 0.0003 

± 0.0002 


d 

0.61682 

20.3653 

0.821 


± 0.00007 

± 0.0003 

± 0.0002 

S2 

X 2 

2.6872 

1.6023 


± 0.0008 

± 0.0002 


1.2091 

1.4031 


± 0.0004 

± 0.0001 


0 

1.2887 


± 0.0001 


d 

1.429 


± 0.0001 


magnetic parameters 


S3 

X3 

h 

S 3 

~ 0 

0.29928 

7.27931 

0.02444 


±0.00005 

±0.00005 

±0.00009 


1.2091 

0.2937 

7.18931 

0.09179 


±0.0004 

±0.00005 

±0.00005 

±0.00008 


0 

0.31041 

7.20405 

0.0712 


±0.00005 

±0.00005 

±0.00009 


d 

0.26302 

7.15342 

0.11427 


±0.00005 

±0.00005 

±0.00009 

S4 

X4 

~ 0 

-0.30974 


±0.00006 


1.2091 

-0.26398 


±0.0004 

±0.00004 


0 

-0.25707 


±0.00004 


d 

-0.42512 


±0.00004 


seems to be a non-physical result, but perhaps a gener¬ 
ous interpretation of effective polarization would allow 
this. The second problem is more serious. Away from 
the unit-cell center, where only the constant magnetic 
susceptibility slab is present, the wave speed exceeds the 


speed of light in vacuum, c. This is the only case where 
this violation is found. In all other models, and all other 
unit-cells analyzed, the frequency independent response 
supports only evanescent waves or waves with speeds less 
than c, or has zero thickness. When this is the case, the 
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FIG. 5. 2D isotropic negative index unit cell. From simulation, the real and imaginary part of the (a) reflection coefficient and 
(b) transmission coefficient, (c) Best global fits of the four models to the simulated data (black dots). The narrow range of 
frequencies displayed is indicated by the dashed box in (a). The models are shown in: green (multi-thickness), yellow (single 
thickness), red (thin sheet) and blue (homogeneous), (d) The combined residuals, A, as defined in equation (261, with the 
same color scheme. 


TABLE II. 2D isotropic negative index unit cell. 


multi-thickness 

xVdof 1 

probability 1 

parameters 14 

single thickness 

500 
~ 0 

11 

thin sheet 

700 
~ 0 

11 


homogeneous 

700 
~ 0 

8 

electric parameters 

Si 

0 


1.4122 

±0.0004 

0 




d 


X? 

0.321 

±0.0006 

0.4367 

±0.0001 

0.2938 


±0.00008 


0.5606 

±0.0001 

/i 

12.8644 

±0.0004 

13.2675 

±0.00009 

13.3235 


±0.00008 


13.2372 

±0.00009 


0.065 

±0.0001 

0.06952 

±0.00007 

0.10551 


±0.00009 


0.05856 

±0.00007 

S2 

0 


1.4122 

±0.0004 

0 




— 


X2 

0.7124 

±0.0007 

0.342 

±0.0003 

0.748 


±0.0004 


— 


/2 

10.8788 

±0.0001 

9.896 

±0.0006 

9.9264 


±0.0007 


— 


^2 

0.0547 

±0.0001 

1.315 

±0.001 

3.344 


±0.003 


- 


S3 

2.2578 

±0.0005 

1.4122 

±0.0004 

0 




d 


X3 

7.4204 

±0.0007 

6.9715 

±0.0006 

5.7128 


±0.0004 


8.2337 

±0.0005 

magnetic parameters 

S4 

0.264 

±0.002 

1.4122 

±0.0004 

0 




d 


X4 

0.3748 

±0.0001 

0.32425 

±0.00004 

0.40891 


±0.00004 


0.26485 

±0.00003 

U 

11.9281 

±0.0009 

10.8916 

±0.00005 

10.9127 


±0.00005 


10.9148 

±0.00005 

54 

0.0691 

±0.0001 

0.06142 

±0.00006 

0.02399 


±0.00006 


0.14302 

±0.00007 

S5 

0.044 

±0.004 

1.4122 

±0.0004 

0 




d 


Xs 

-0.4018 

±0.0003 

-0.3861 

±0.00007 

-0.5123 


±0.0001 


-0.47033 

±0.00005 











































frequency independent response is causal, and it remains 
so when combining with a causal Lorentzian response. 
It may be that the fit parameters could be constrained 
to prevent causality violations, and compelling fits still 
obtained, though that was not confirmed here. 


IV. STATISTICAL MODEL SELECTION 


To compare different models for a given data set, one 
must quantify the goodness of fit as well as assess a 
penalty for each free fit parameter. This penalty is some¬ 
times called the Occam’s factor. An expression for the 
posterior probability of model correctness can be found 
using Bayes Theorem. We follow the notation and con¬ 
ventions of the book by Sivia and Skillin^Sl. 

prob(J|M)prab(M) 

prob [Dj 

This expression gives the posterior probability that the 
model, M, is correct given the data, D. The required fac¬ 
tors include: the probability of the data given the model 
(also called the likelihood function), the prior probabil¬ 
ity for the model, and the probability of the data. The 
last is usually not directly calculable, since it would re¬ 
quire integrating over all possible models. The need for 
it can be eliminated by using a normalization condition, 
or by seeking only a ratio for model-correctness probabil¬ 
ities between two alternative models. Here the data, D, 
refers to the reflection and transmission coefficients from 
a simulation. If we assume that the log of the parameter 
dependent likelihood function can be well approximated 
by a quadratic series expansion around the best fit pa¬ 
rameter values, Aq, and use a uniform probability over a 
finite interval for the prior parameter values, we find the 
parameter dependent likelihood function to be, 

prob {D\M) = prob (£>1 Aq, M) (19) 

i 

where the 5Xi are the parameter uncertainties and the 
AAi are the prior parameter ranges. The product of fac¬ 
tors on the right comprise the Occam’s penalty for adding 
free parameters. (The series approximation is the same 
used in the standard calculation of the covariance ma¬ 
trix.) The ratio of model correctness probabilities for 
two models. Mi and M 2 is then 


prob (Ml ID) 
prob (M 2 \D) 


prob(D| Ao,Mi)n'^ 

i 

prob {D I /xo, M 2 ) n 

3 ^ 


( 20 ) 


where the Xi are the parameters for model Mi and the 
fjLj are the parameters for model M 2 . This assumes that 
the two models are equally likely prior to examining the 
data 


prob (Ml) ^ ^ 
prob (M 2 ) 


( 21 ) 


The parameter uncertainties are given by the square-root 
of the diagonal elements of the covariance matrix. 

5X, = ( 22 ) 

where as usual, the covariance matrix is given by the 
inverse of the Hessian (matrix of second derivatives) of 
the log of the likelihood function evaluated at the best 
fit parameters. 

S = -(VaVaT(Ao))-^ (23) 

For a least-squares fit, the log of the likelihood function 
is 

L (A) = constant — (A) (24) 

where is the sum of the squares of the normalized 
residuals 


(A) 





(25) 


where 

A.(A) = 7|rf(A)-rfe|V|t'^(A)-t'|^ (26) 

fc is a frequency index, r^(A) and t'^^ {\) comprise the 
model evaluated at frequency index k and fit parameters 
A, and and t'j. are the data. We have assumed that 
the simulation data uncertainty, ct, is independent of fre¬ 
quency, and the same for all of the real and imaginary 
parts of Tfc and t).. Note that in the parameter dependent 
likelihood probability 

prob(D I Ao,M) = exp(L(Ao)) oc exp ^-^X^ (Ao) 

(27) 

the proportionality constant does not depend on the 
model and may be neglected when computing the model 
probability ratios. 

There are several issues that arise with the data un¬ 
certainty, tJ, when the data is generated in a simulation. 
The least-squares minimization procedure is derived un¬ 
der the assumption that cr is known, and describes the 
width of an independent Gaussian stochastic variable. 
However, the error present in FDTD simulation data is 
all systematic and is usually dominated by finite mesh 
effects and oscillations introduced in the frequency do¬ 
main variables by truncation of the transient variable re¬ 
sponse. We take the viewpoint that a simulation with a 
given mesh is a linear system with a valid response func¬ 
tion. Thus, models may be compared using simulation 
data with course or fine meshes. (Of course, to best ap¬ 
proximate a “continuum” physical system, convergence 
of the response with respect to mesh density should be 
sought.) The remaining source of error, the transient- 
truncation induced oscillation, is, unfortunately, neither 
independent nor stochastic. 
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FIG. 6. ELC unit cell. From simulation, the real and imaginary part of the (a) reflection coefficient and (b) transmission 
coefficient, (c) Best global fits of the four models to the simulated data (black dots). The narrow range of frequencies displayed 
is indicated by the dashed box in (a). The models are shown in: green (multi-thickness), yellow (single thickness), red (thin 
sheet) and blue (homogeneous), (d) The combined residuals, A, as defined in equation (|26|l, with the same color scheme. 


TABLE III. ELC unit cell. 



multi-thickness 

single thickness 

thin sheet 


homogeneous 

xVdof 

1 

4 

40 


100 

probability 

1 

~ 0 

~ 0 


~ 0 

parameters 

6 

6 

5 


5 


electric parameters 


Si 

X? 

h 

Si 

0 

0.3909 

15.3043 

0.2651 


± 0.0001 

± 0.0003 

± 0.0004 


1.7125 

0.4316 

15.156 

0.2767 


± 0.0009 

± 0.0001 

± 0.0003 

± 0.0003 


0 

0.33828 

15.2831 

0.3432 


± 0.0001 

± 0.0003 

± 0.0004 


d 

0.5388 

15.0181 

0.4076 


± 0.0002 

± 0.0003 

± 0.0003 

S2 

X2 

0 

1.872 


± 0.0005 


1.7125 

2.279 


± 0.0009 

± 0.0004 


0 

1.798 


± 0.0002 


d 

2.8887 


± 0.0004 


magnetic parameters 


S3 

5.114 

± 0.003 

1.7125 

± 0.0009 

0 


d 


X3 

- 0.8308 

± 0.0009 

- 0.25693 

± 0.00008 

- 0.26243 

± 0.0001 

- 0.43796 

± 0.00006 


The lack of independence manifests as a correlation of 
the error over a significant range of frequency. For the 
model probability calculations, we mitigate this problem 
by decimating the data to a courser, evenly-spaced sub¬ 
set. The data are originally very hnely-spaced in fre¬ 
quency, for accurate model fitting. By so decimating. 
We were able to reduce the (off-peak) auto-correlation of 
the residuals by a factor of five, while still capturing the 


significant behavior of the reflection and transmission co¬ 
efficients. For the SRR, 2D isotropic and ELC unit-cells, 
the original number of frequency points was reduced from 
10000, 10000, and 12500 to 121, 130, and 39 respectively. 

The lack of stochasticity means that we cannot es¬ 
timate the magnitude of the error by making multiple 
measurements (i.e. simulations). However, due the re¬ 
markable quality of the fits using the multiple-thickness 
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model, we can take this model to be a proxy for the un¬ 
known exact results, and use the residuals of the fits to 
this model to estimate a. With this estimate, we can 
test our assumption that a is independent of frequency, 
and the same for all of the real and imaginary parts of 
rfc and We find that real and imaginary parts of rk 
and do indeed have very similar residuals, but the fre¬ 
quency independence is more questionable. Our method 
of estimating a normalizes the / DOF and model prob¬ 
ability to unity for the multiple-thickness model, as seen 
in Tables E0 and [ml The four models can then only be 
judged in a relative sense. 

Finally, a slight complication arises when the off- 
diagonal elements of the covariance matrix are not negli¬ 
gible, which can result in an underestimation of the pa¬ 
rameter uncertainties. One can correct this by finding 
a new set of parameters that diagonalize the covariance. 
Since the covariance is a real symmetric matrix, the di¬ 
agonalizing parameters are given by 


A' = U^A (28) 

where the normalized Eigenvectors of the covariance com¬ 
prise the columns of the orthogonal matrix U. The new 
parameter uncertainties are given by 

^A' = ^/^■ (29) 


by 




(30) 


For the likelihood function, (19), we had assumed that 


the a-priori range probabilities were uniform. This is no 
longer the case, but we neglect correcting for this for 
simplicity. The model probabilities in Tables El and 
|III| are calculated using the diagonalized forms, but the 
parameter values and parameter uncertainties refer to the 
original, physically-derived model parameters. 


V. CONCLUSION 

We have presented four models for representing the 
response of metamaterials to normally incident plane 
waves. All of these models rely on causally-responding 
components (with a, perhaps correctable, exception 
noted above). Fitted models can thus provide a com¬ 
parison of unit-cell designs, using physically meaningful 
figures of merit. Only one of the models—the multi¬ 
thickness model—exhibits compelling, and quantitatively 
precise representation of the simulated reflection and 
transmission behavior for the three unit-cells here an¬ 
alyzed. We believe this to be the simplest model thus far 
presented in the literature that meets the criteria of pro¬ 
viding physically meaningful figures of merit, and quan¬ 
titatively precise representation. 
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